# -*- coding: utf-8 -*-
import numpy as np
import pandas as pd

class AerodynamicCalculator:
    def __init__(self, U, a, rho=1.225):
        """气动参数计算器
        参数：
            U: 来流速度 (m/s)
            a: 圆柱半径 (m)
            rho: 空气密度 (kg/m³)，默认1.225
        """
        self.U = U
        self.a = a
        self.rho = rho

    def compute_pressure(self, z_samples):
        """计算压力系数分布
        参数：
            z_samples: 复数坐标数组 (np.ndarray[complex])
        返回：
            Cp: 压力系数数组 (np.ndarray)
        """
        V = self.U * (1 - (self.a**2) / (z_samples**2))
        V_mag = np.abs(V)
        Cp = 1 - (V_mag / self.U)**2
        return Cp

    def export_report(self, filename="pressure_report.csv"):
        """生成标准化气动分析报告
        参数：
            filename: 输出文件名 (str)
        """
        theta = np.linspace(0, 2*np.pi, 36)  # 10°间隔采样
        z_samples = self.a * np.exp(1j * theta)
        Cp = self.compute_pressure(z_samples)
        V = self.U * (1 - (self.a**2)/(z_samples**2))

        # 构建DataFrame
        df = pd.DataFrame({
            'theta_deg': np.degrees(theta),
            'Cp': Cp,
            'V_real': V.real,
            'V_imag': V.imag
        })
        
        # 数据保存
        df.to_csv(filename, index=False, float_format='%.4f')
        print(f"报告已保存至 {filename}")

    def validate_pressure(self, tolerance=1e-3):
        """压力系数物理验证
        验证驻点(0°和180°)和极值点(90°和270°)的压力系数
        """
        test_angles = np.array([0, np.pi/2, np.pi, 3*np.pi/2])
        z_test = self.a * np.exp(1j * test_angles)
        Cp_test = self.compute_pressure(z_test)
        
        # 理论值验证
        assert np.allclose(Cp_test[[0, 2]], [1, 1], atol=tolerance), "驻点压力验证失败"
        assert np.allclose(Cp_test[[1, 3]], [-3, -3], atol=tolerance), "极值点压力验证失败"
        print("所有物理验证通过")

# 使用示例
if __name__ == "__main__":
    # 初始化计算器（参数与GUI默认值一致）
    calc = AerodynamicCalculator(U=5.0, a=1.0)
    
    # 执行物理验证
    calc.validate_pressure()
    
    # 生成报告
    calc.export_report()
    
    # 压力分布可视化（可选）
    import matplotlib.pyplot as plt
    theta = np.linspace(0, 2*np.pi, 36)
    Cp = calc.compute_pressure(calc.a * np.exp(1j * theta))
    
    plt.figure(figsize=(8, 8))
    ax = plt.subplot(111, polar=True)
    ax.plot(theta, Cp, 'r-', lw=2, label='计算结果')
    ax.plot(theta, 1 - 4*np.sin(theta)**2, 'b--', lw=1, label='理论解')  # 理论曲线
    ax.set_title(f'圆柱表面压力分布 (U={calc.U}m/s, a={calc.a}m)')
    ax.legend(loc='best')
    plt.savefig('pressure_distribution.png')
    plt.close()
    
    

